home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / fft1 / fast.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  6.6 KB  |  324 lines
  1. /* *****************************************************************************
  2. *
  3. * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  4. * All Rights Reserved.
  5. *
  6. * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  7. * the contents of this file may not be disclosed to third parties, copied or
  8. * duplicated in any form, in whole or in part, without the prior written
  9. * permission of Silicon Graphics, Inc.
  10. *
  11. * RESTRICTED RIGHTS LEGEND:
  12. * Use, duplication or disclosure by the Government is subject to restrictions
  13. * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  14. * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  15. * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  16. * rights reserved under the Copyright Laws of the United States.
  17. *
  18. ***************************************************************************** */
  19. #include <gl.h>
  20. #include <device.h>
  21. #include <math.h>
  22. #include "fast.h"
  23. #include "event.h"
  24. #include "fft.h"
  25.  
  26. void display();
  27. void InitFast();
  28. void matgen();
  29. void dir_fft(), inv_fft();
  30. double second();
  31.  
  32.  
  33. extern int mainmenu, realmenu, fouriermenu, fourier2menu, quitmenu;
  34.  
  35. int do_exit;
  36. int ask_quit;
  37.  
  38. int x_size, y_size;
  39. int this_proc;
  40.  
  41. float *mat, *save, *power;
  42. int * mycolors;
  43. int *ipvt, lda, info, x_0, y_0;
  44.  
  45. double this_time, this_flop;
  46.  
  47. char title[128];
  48.  
  49.  
  50. fast (progname)
  51. char *progname;
  52. {
  53.     InitFast();
  54.     Init(progname);
  55.  
  56.     do_exit = FALSE;
  57.  
  58.     while(! do_exit) {
  59.         event();
  60.     }
  61. }
  62.  
  63. do_frame(){
  64.     
  65.     char str[32];
  66.     float v[2];
  67.  
  68.    viewport(0,WIN_SIZE+WIN_RIGHT,0,WIN_SIZE);
  69.  
  70.     cpack (oBLACK);
  71.     clear ();
  72.  
  73.      cpack (oGREEN);
  74.     ortho2(0,WIN_SIZE+WIN_RIGHT,0, WIN_SIZE);
  75.  
  76.     sprintf( str, "SIZE = %d x %d", x_size, y_size);
  77.     cmov2i( 250, WIN_SIZE - 30);
  78.     charstr(str);
  79.  
  80.     cmov2i( 100, 20);
  81.     charstr(title);
  82.  
  83.     bgnpolygon();
  84.     v[0] = WIN_SIZE;
  85.     v[1] = WIN_SIZE/2;
  86.     cpack(oBLACK);
  87.     v2f(v);
  88.     v[1] -= 100;
  89.     cpack(0x007fff);
  90.     v2f(v);
  91.     v[0] += 50;
  92.     v2f(v);
  93.     v[1] += 100;
  94.     cpack(oBLACK);
  95.     v2f(v);
  96.     endpolygon();
  97.     v[1] -= 100;
  98.     bgnpolygon();
  99.     cpack(0x007fff);
  100.     v2f(v);
  101.     v[1] -= 100;
  102.     cpack(0x00ffff);
  103.     v2f(v);
  104.     v[0] -= 50;
  105.     v2f(v);
  106.     v[1] += 100;
  107.     cpack(0x007fff);
  108.     v2f(v);
  109.     endpolygon();
  110.  
  111.     bgnpolygon();
  112.     v[0] = WIN_SIZE;
  113.     v[1] = WIN_SIZE/2;
  114.     cpack(oBLACK);
  115.     v2f(v);
  116.     v[1] += 100;
  117.     cpack(0xff7f00);
  118.     v2f(v);
  119.     v[0] += 50;
  120.     v2f(v);
  121.     cpack(0xff7f00);
  122.     v2f(v);
  123.     v[1] -= 100;
  124.     cpack(oBLACK);
  125.     v2f(v);
  126.     endpolygon();
  127.     v[1] += 100;
  128.     bgnpolygon();
  129.     cpack(0xff7f00);
  130.     v2f(v);
  131.     v[1] += 100;
  132.     cpack(oCYAN);
  133.     v2f(v);
  134.     v[0] -= 50;
  135.     v2f(v);
  136.     v[1] -= 100;
  137.     cpack(0xff7f00);
  138.     v2f(v);
  139.     endpolygon();
  140.     cpack(oWHITE);
  141.     rect(WIN_SIZE,WIN_SIZE/2-200,WIN_SIZE+50,WIN_SIZE/2+200);
  142.  
  143.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 - 200 -4);
  144.     charstr(" -Max");
  145.  
  146.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 + 200 -4);
  147.     charstr(" +Max");
  148.  
  149.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 -4);
  150.     charstr(" 0.0 ");
  151.  
  152. }
  153. ShowFast()
  154. {
  155.     x_0 = (WIN_SIZE/2)-(x_size/2);
  156.     y_0 = (WIN_SIZE/2)-(y_size/2);
  157.     do_frame();
  158.     viewport( x_0, x_0 + x_size + 2, y_0, y_0 + y_size + 2);
  159.     ortho2(-1, x_size+1,-1, y_size+1);
  160.     DrawFast();
  161. }
  162.  
  163. void do_quit()
  164. {
  165.     if(mat)
  166.         free(mat);
  167.     gexit();
  168.     exit(0);
  169. }
  170.  
  171.  
  172. void do_rightmouse()
  173. {
  174.     ask_quit = 0;
  175.  
  176.     dopup(mainmenu);
  177.  
  178.     if(ask_quit)
  179.     dopup(quitmenu);
  180.  
  181. }
  182.  
  183. void display()
  184. {
  185.     do_frame();
  186.     swapbuffers();
  187.     ShowFast();
  188.     swapbuffers();
  189. }
  190.  
  191. void do_resize(){
  192.     
  193.     lda = y_size + 1;
  194.  
  195.     if( mat){ 
  196.     mat  = (float *)realloc( mat, lda * x_size * sizeof(float));
  197.     power  = (float *)realloc( power, lda * x_size * sizeof(float));
  198.     mycolors  = (int *)realloc( mycolors, lda * x_size * sizeof(int));
  199.     save  = (float *)realloc( save, (y_size + 3*x_size + 4 * FACTOR_SPACE) * sizeof(float));
  200.     }
  201.     else{ 
  202.     mat   = (float *)malloc( lda * x_size * sizeof(float));
  203.     power   = (float *)malloc( lda * x_size * sizeof(float));
  204.     mycolors   = (int *)malloc( lda * x_size * sizeof(int));
  205.     save  = (float *)malloc( (y_size + 3*x_size + 4 * FACTOR_SPACE) * sizeof(float));
  206.     }
  207.  
  208.     sfft2di( y_size, x_size, save);
  209.     matgen();
  210.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  211.     mainmenu = realmenu;
  212.     strcpy( title, "XY domain");
  213.     display();
  214. }
  215. void add_noise() {
  216.  
  217.     matnoise_( mat, &lda, &y_size, &x_size);
  218.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  219.     sprintf(title,(mainmenu == realmenu) ? "XY domain" : "FOURIER domain");
  220.     display();
  221. }
  222.  
  223. void this_pow(){
  224.     matpower_(&y_size, &x_size, mat, power, &lda);
  225.     matlog_(&y_size, &x_size, power, &lda);
  226.     matcolor_(&y_size, &x_size, power, mycolors, &lda);
  227.     sprintf(title,"Power Spectrum (Log Scale)");
  228.     display();
  229. }
  230.  
  231. void size500_500(){
  232.  
  233.     x_size = 500;
  234.     y_size = 500;
  235.     do_resize();
  236. }
  237.  
  238. void size512_512(){
  239.  
  240.     x_size = 512;
  241.     y_size = 512;
  242.     do_resize();
  243. }
  244.  
  245. void size576_405(){
  246.  
  247.     x_size = 576;
  248.     y_size = 405;
  249.     do_resize();
  250. }
  251.  
  252. void dir_fft(){
  253.  
  254.     this_time = second();
  255.     sfft2d( -1, y_size, x_size, mat, lda, save);
  256.     this_time = second() - this_time;
  257.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  258.     mainmenu = fouriermenu;
  259.     sprintf(title,"Direct (%d x %d) FFT performed in %.2f Seconds \n",
  260.         y_size,x_size, this_time);
  261.     display();
  262. }
  263. void inv_fft(){
  264.     float alpha = 1./(float)(x_size*y_size);
  265.  
  266.     this_time = second();
  267.     sfft2d( 1, y_size, x_size, mat, lda, save);
  268.     sscal2d_( &y_size, &x_size, &alpha, mat, &lda);
  269.     this_time = second() - this_time;
  270.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  271.     mainmenu = realmenu;
  272.     sprintf(title,"Inverse (%d x %d) FFT performed in %.2f Seconds \n",
  273.         y_size,x_size, this_time);
  274.     display();
  275. }
  276.  
  277. void this_quit(){
  278.  
  279.     ask_quit = 1;
  280. }
  281. void InitFast(){
  282.  
  283.     this_proc = mp_numthreads_();
  284.     x_size = 512;
  285.     y_size = 512;
  286. }
  287. #define SCALE        20
  288. #define XLENGTH        19*SCALE
  289. #define YLENGTH        11*SCALE
  290.  
  291. void myrec( float *pt, int xl, int yl, double val){
  292.   int i, j;
  293.     for( j = 0 ; j < xl; j++, pt += lda)
  294.     for( i = 0; i < yl ; i++)
  295.         pt[i] += val;
  296. }
  297. myF( float *pt, double val){
  298.     myrec( pt, SCALE, 7*SCALE, val);
  299.     myrec( pt+SCALE*lda+3*SCALE, 2*SCALE, SCALE, val);
  300.     myrec( pt+SCALE*lda+6*SCALE, 3*SCALE, SCALE, val);
  301. }
  302. myT( float *pt, double val){
  303.     myrec( pt+0*SCALE*lda+6*SCALE, 5*SCALE, 1*SCALE, val);
  304.     myrec( pt+2*SCALE*lda+0*SCALE, 1*SCALE, 6*SCALE, val);
  305. }
  306. void mywrite( float *pt){
  307.     myrec( pt, XLENGTH, YLENGTH, 1.);
  308.     myF( pt + 2 * SCALE +  2 * SCALE*lda, -2.);
  309.     myF( pt + 2 * SCALE +  7 * SCALE*lda, -2.);
  310.     myT( pt + 2 * SCALE + 12 * SCALE*lda, -2.);
  311. }
  312.  
  313. void matgen(){ 
  314.     bzero( mat,  lda * x_size * sizeof(float));
  315. }
  316. void add_box(){
  317.     int xor, yor;
  318.     xor = x_size/2 - XLENGTH/2;
  319.     yor = y_size/2 - YLENGTH/2;
  320.     mywrite( mat+yor+xor*lda);
  321.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  322.     display();
  323. }
  324.